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The increasing complexity of heat 
transfer and process situations which 
involve fluid flow has demanded the 
frequent use of flow passages of unusual 
geometrical configuration. The present 
investigation is concerned with one such 
novel configuration, namely the longi- 
tudinal flow between solid cylindrical 
rods which are arranged in regular array. 
A schematic diagram of the situation 
under study is shown in Figure 1. The 
rods may be located either in triangular 
or square array. The flow will be taken 
to be laminar and fully developed. 

The aim of this analysis is to determine 
the pressure drop, shear stress, and 
velocity-distribution characteristics of 
the system. The starting point of this 
study is the basic law of momentum 
conservation. The resulting differential 
equation has been solved in an approxi- 
mate, but almost exact, manner by the 
use of truncated trigonometric series. 
Results are obtained over a wide range of 
porosity values for both the triangular 
and square arrays. Heat transfer has not 
been considered. 

The configuration under investigation 
has potential application in compact heat 
exchangers for nuclear reactors and other 
situations. Further the results should 
also be of interest in the theory of flow 
through unconsolidated porous beds (la, 
2a). 

The only related analytical work 
known to the authors is that of Emersle- 
ben (5), who considered only the square 
array. His rather involved solution, based 
on complex zeta functions, appears to be 
valid only at high porosities. Experiments 
covering a porosity range of 0.093 to 
0.984 have been made by Sullivan (4) 
using parallel-oriented fibers, most of 


the tests being for fibers in random array. 
These previous investigations will be 
compared with the present theory in a 
later section. 

ANALYSIS 

The Governing Equation and 
Its General Solution 

The physical principle which governs 
the velocity distribution in a (isothermal) 
flowing fluid is conservation of momen- 
tum. To translate this physical law into 
mathematical terms, cylindrical coordi- 
nates will be used. The derivation is 
facilitated by Figure 2. The conser- 
vation principle requires that under 
steady state conditions the net change of 
momentum must equal the net forces. 
However for fully developed flow the net 
momentum flux is zero (no accelerations), 
and hence the forces must sum to zero. 
The forces involved in the problem are 
those of pressure and viscosity, and these 
must be in balance. Therefore 


ential equation is frequently called 
Poisson's equation. 

It is not difficult to write a solution of 
this equation, but subsequently one must 
face the more challenging task of fitting 
the solution to the particular boundary 
conditions of specific problems. The 
general solution is dealt with here, leaving 
specific flow situations to later sections. 

In approaching Equation (2) it is 
convenient to introduce a reduced ve- 
locity defined by 



By substitution into Equation (2) one 
finds that u * mus£ obey 


aV . l l aV _ 
dr 2 + r dr + r 2 d0 2 


( 4 ) 


which is the well-known Laplace equa- 
tion. The general solution of Equation 
(4) may be taken from numerous books 
on advanced calculus. Therefore 


[(p + g&)- p } d$ dr 

-»[%(%' dtd ‘) dr (i > 
+ fi{}fe drd ‘) d<l ] 

When one rearranges, the governing 
equation for the velocity as a function of 
the coordinates r and 6 is obtained: 

d 2 u . 1 du . 1 d 2 u 1 dp ( . 
dr 2 ^rdr^r 2 dd 2 "pdz W 

where the pressure gradient dp/dz is a 
negative constant. This partial differ- 


u* = A + B In r 

+ £ (<V + D k r~ k ) ( 5 ) 

*-t 

•(. E k cos kd + F k sin kd) 

where k takes on integral values to 
ensure that the velocity is single valued, 
that is that the velocity computed at a 
location (r, 0) is identical to that com- 
puted at (r, 6 + 2t). The constants A , 
B, etc., in Equation (5) are to be de- 
termined from the boundary conditions, 
as are the number of terms of the series. 

Thus a general solution for the velocity 
is obtained by the combining of Equa- 
tions (3) and (5) which yields 
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Fig. 1. Schematic of flow configuration. 



Fig. 2. Control volume for deriving momen- 
tum equation. 


u=A+B lar-'j(-if) 

+ Z (C k r k + D k r ~ k ) (6) 

k~\ 

•( E k cos kQ + F k sin kO) 

Now one may turn to the problem of 
specializing this solution to flow parallel 
to cylindrical rods in regular array. 

Cylinders in Triangular Array 

Consideration is first given to flow 
between cylinders arranged on centers 
which form the vertices of equilateral 
triangles. An end view of the arrangement 
is shown in Figure 3. From the symmetry 
of the situation it is easily seen that 
attention need be focused on only the 
cross-hatched element in the left-hand 
sketch. An enlarged view is shown at the 
right, on which has been noted the 
boundary conditions and dimensional 
nomenclature. The condition du/dn = 0 
is the expression of the symmetry 
property, while the condition that u = 0 


Table 1. Listings of A, and 5, foe Various Spacing Ratios 


(a) Values of A,-. 


s/r„ 

Ai 

A2 

A* 

A* 

A* 

A* 

A 7 

4.0 

- 0.0505 

- 0.0008 

0.0000 




► 

2.0 

- 0.0505 

- 0.0008 

0.0000 




► 

1.5 

- 0.0502 

- 0.0007 

0.0000 




► 

1.2 

- 0.0469 

0.0007 

0.0002 

0.0000 



> 

1.1 

- 0.0416 

0.0028 

0.0004 

0.0000 



> 

1.05 

- 0.0368 

0.0043 

0.0003 

,- 0.0001 

0.0000 


> 

1.04 

- 0.0357 

0.0046 

0.0002 

- 0.0001 

0.0000 


> 

1.03 

- 0.0345 

0.0049 

0.0002 

- 0.0001 

0.0000 


> 

1.02 

- 0.0332 

0.0051 

0.0000 

- 0.0001 

0.0000 


> 

1 . 01 - 

- 0.0319 

0.0052 

- 0.0001 

- 0.0002 

0.0000 


> 

1.00 

- 0.0305 

0.0053 

- 0.0003 

- 0.0002 

0.0000 ■ 


— 





( 6 ) Values of 5 , 




s/ro 

5i 


$3 



5 « 


4.0 

- 0.1253 

- 0.0106 

- 0.0006 

0.0000 



— 

2.0 

- 0.1250 

- 0.0105 

- 0.0006 

0.0000 



> 

1.5 

- 0.1225 

- 0.0091 

- 0.0002 

0.0000 



> 

1.2 

- 0.1104 

- 0.0024 

- 0.0015 

0.0003 

0.0001 

0.0000 

> 

1.1 

- 0.0987 

0.0036 

0.0029 

0.0005 

0.0000 


> 

1.05 

- 0.0904 

0.0073 

0.0032 

0.0002 

- 0.0001 

0.0000 

> 


at the inner boundary stems from the 
no-slip requirement of viscous flow. 

One may now proceed to apply the 
boundary conditions to determine the 
constants of Equation (6). Starting with 
the simplest conditions, one requires 
first that du/dn = du/dO = 0 at 6 = 0 
deg. and at 0 = 30 deg. From the first 
of these it follows that since cos 0^0, 
then 

. F k ee 0 (7a) 

while from the second one finds that 


k = 6, 12, 18- •• (76) 

to guarantee that sin for/6 = 0. Next 
imposing the condition that u = 0 at 
r == r 0 one gets 


D k = -<V, 2 \ 

A = —B In r 0 + 


r 0 


(_l&) 

\ n dz) 


(7c) 


Further it is required that the total drag 
force exerted on the fluid by the solid rod 
be balanced by the net pressure force 
acting over the entire cross section of the 
typical element (Figure 36); that is 



The evaluation of this over-all force 
balance from Equation (6) yields 



Before going on to the final boundary 
condition the findings of the previous 
paragraph are brought together, and 
Equation (6) becomes 



where G , = C/E,-. There still remains the 
task of determining the G iy and at one’s 
disposal is the condition that du/dn = 0 
on the right-hand boundary of Figure 
36, on which r = s / cos 6. It is convenient 
to make use of the identity 


du 

dx 


du 

— ~ cos 6 — 
dr 


du sin 0 
dd r 


where one may associate x with n on the 
boundary under consideration. Intro- 
ducing Equation (9) and setting 
du/dn = 0 when r = s / cos 6 , one finds 
after rearrangement 


X A,- (cos 0) 1 4 cos (6 j — 1)0 


+ (^)'”co 5 < | y + l) l >] (7e> 



Equation (7e) provides a means for 
determining the A, (that is (?,-). The first 
thought for attacking this equation 
would be to apply the techniques of 
Fourier series. Such an approach, if 
possible, would provide an infinite set 
of A,- while satisfying Equation (7e) at 
all points along the boundary (that is all 
O°<0< 30°). Unfortunately the nature 
of Equation (7e) precludes the use of 
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(a) end vie* 


(b) typical eleoent 


Fig. 3. Diagram of equilateral triangular array. 


Fourier analysis, and some other ap- 
proach must be found for solving for the 

A,, 

The method used here is to apply 
Equations (7e) at a finite number of 
points along the boundary. Suppose for 
example that one successively evaluates 
Equation ( 7e ) at six values of 0 between 
0 and 30°. This will provide six equations. 
If one truncates the series after six terms 
(that is l — 6), then there are sufficient 
simultaneous equations to evaluate the 
coefficients Ai, A 2 • • • A®. This procedure 
may now be repeated with the use of 
seven boundary points and seven series 
coefficients, and so forth. The sets of 
coefficients A, from these repeated calcu- 
lations were compared. It was im- 
mediately seen that adding additional 
terms to the series did not significantly 
affect the numerical values of the first 
few coefficients. Further it was found 
that only these first coefficients are 
important in the computation of the 
shear stress and velocity distribution. 
In view of these favorable circumstances 
the method presented here for determin- 
ing the A, appears quite adequate. 

Numerical values of A,-, computed as 
outlined above, have been listed in Table 
1. Inasmuch as the ratio s/r 0 (half 
spacing to rod radius) appears in Equa- 
tions (7e), so does it appear as a param- 
eter of Table 1. It is noted that for large 
spacings (that is large s/r 0 ) the A,- are 
little affected by increases of spacing, as 
might have been expected since s/r 0 is 
raised to a large negative power in 
Equation (7e). The tabulation is given 
to four decimal places because this is 
sufficient for the shear stress and velocity 
computations.* 

So with the determination of the A,- 
(that is, Gi ) one may return to Equa- 
tion (9) and state that the velocity dis- 
tribution for the triangular array is now 
available. In a later section the authors 
will make use of this velocity solution to 
compute several quantities of engineering 
interest. But first they will solve for the 
velocity distribution associated with flow 
between rods in square array. 

Cylinders in Square Array 

An end view of the configuration for 
flow between cylinders in square array is 

♦Additional figures were used in satisfying Equa- 
tion (7c). 


shown schematically in Figure 4. Because 
of symmetry considerations attention is 
confined to the cross-hatched element of 
the sketch. This element has essentially 
the same form as that shown in Figure 
36, except that now the opening angle is 
45 deg. rather than 30 deg. as before. 

The velocity solution as given by 
Equations (6) applies, and one now has 
to find the constants A, £, • * * appro-, 
priate to the square array. The pro- 
cedure is identical to that followed in 

the preceding section for the triangular 
array. From the conditions that du/dn — 
du/dd = 0 attf = 0 deg. and 0 = 45 deg., 
there is obtained, respectively 

F k = 0 (10a) 

and 

* = 4 , 8 , 12 — - ( 106 ) 

Imposing the requirement that u = 0 
at r = r 0 gives 


D k = — C*r 0 2 *, 


A = 


-B In r„ + ^ ( 



(10c) 


while the over-all force balance between 
net pressure and wall shear, Equation 
(8) with tt/ 6 replaced by 7 t/ 4, provides 
the following value of B : 



Introducing these findings into (6) 
results in the following equation for the 
velocity: 



where the G , (= C, •!£,■) still remain to be 
determined. The condition that du/dn = 
0 on the right-hand boundary of Figure 
36, on which r = s / cos 6 , may be used. 
The procedure for determining the (?, 
follows along the lines previously dis- 
cussed. For several values of 6 between 






Fig. 4. Diagram of square array. 



Fig. 5. Diagram for pressure-drop and 
friction-factor derivation. 


0 and 45 deg. one successively applies 
the condition of zero normal derivative; 
that is 


i 

£ 6/ (cos 0) 1 ' 4 ' cos (4; — 1)0 
j-i L 


+ 


^ r 0 cos 0 y* 


cos (4 j + 1) 0 




(12a) 



The result is l simultaneous equations for 
6i, 02 * * * 0i- Sufficient terms are retained 
in the series to assure good accuracy in 
shear-stress and velocity-distribution 
calculations. Numerical values of the 0,- 
obtained in this manner are listed in 
Table 1 as a function of s/r 0 . 

So with the determination of the 0,* 
(that is G,) the velocity solution for the 
square array, Equation (11), can be 
regarded as known, and one can pass on 
to the presentation of results of engineer- 
ing interest. 


RESULTS 

Attention will first be focused on the 
pressure drop-flow relationship, which is 
generally the result of greatest practical 
importance. Then this information will 
be rephrased in terms of the friction 
factor and Reynolds number. Finally 
plots will be given of velocity contours 
and also of the distribution of the wall 
shear stress around the periphery of a rod. 
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Pressure-Drop— Flow Relationship 

First the volume rate of flow Q which 
may be calculated by integrating the 
velocity over the flow area is introduced. 
As before, only a typical element of the 
flow configuration need be used (Figure 
5). In terms of symbols given there the 
integral for Q becomes 


Q -fL 


-n 


ur dr dB 

00 **/C08 0 


(13) 


ur dr dd 


In evaluating Equation (13) one intro- 
duces u from either Equations (9) or (11) 
depending whether a triangular or square 
array is being considered. As is shown in 
the Appendix the end result of inserting 
the velocity expressions takes the form 


Q = 



To 

-^-•function (s/r 0 ) 


(14) 


or 


dz r 0 4 function (s/r 0 ) 

The function of s/r 0 is different for each 
array and has been computed by the 
use of the integrals given in the Appendix. 

Before numerical results are discussed, 
an important fact may be drawn from 
Equation (14a): for a given array with 
fixed-rod diameter and spacing (that is, 
fixed configuration) the pressure drop is 
directly proportional to flow rate. This 
finding is in accord with previous laminar- 
flow analyses. 

As a final prelude to presenting the 
pressure-drop results, the authors chose 
to replace the ratio s/r 0 by the more 
general variable, the porosity, which 
is the fraction of the total cross section 
available to flow. In terms of Figure 5 
one may write 


Ai + A 2 


= 1 


(15) 


$0 

(s/r 0 ) 2 tan 00 


where the facts that ( A x + A z ) is a right 
triangle and Ai is a circular sector have 
been used in deriving the last expression. 
For the triangular array, 0 O = 30 deg., 
and ' 


_ i Vfor 
e 1 6(sAo) 2 (15a) 

while for the square array, 0 O = 45 deg., 
and 


€ 



(156) 


With these expressions, plus the evalu- 
ated integrals of the Appendix, one is 
able to plot the pressure-drop-flow 
relationship as a function of porosity on 
Figure 6. There are several aspects of the 


curves worthy of discussion. First it is 
seen that for the same flow, r 0 and e, the 
triangular array gives rise to a higher 
pressure drop. This might have been 
intuitively expected because each rod in 



Fig. 6. Pressure-drop-flow parameter as a 
function of porosity. 



Fig. 7. Friction factor-Reynolds number as a 
' function of porosity. 




Fig. 8. Representative velocity contour lines 
for equilateral triangular array. 


the triangular array is surrounded by 
more neighbors and hence feels more flow 
disturbance than does a rod in the square 
array. 

It would appear that the two curves of 
the figure have generally the same shape, 
with possibly a growing deviation at 
smaller values of e. Consideration of the 
system suggests that, except for large 
spacings (large e), similar flow perform- 
ance need not necessarily be expected 
from the different arrays. At large 
spacings, where there is little effect of 
neighboring rods, the flow passages of the 
two arrays are almost geometrically 
similar. But at small spacings the flow 
passages of the two arrays are very 
different indeed. For the case of cylinders 
touching ( s/r 0 = 1) the passages of the 
square array are curvilinear squares, and 
those of the triangular array are curvi- 
linear triangles. One would have no 
reason to expect that at small spacings 
the curves for the two arrays should be 
identical in form. (These observations 
will have even greater meaning when the 
friction-factor curves are given.) 

That the curves of Figure 6 should give 
increasing ordinate values with decreas- 
ing porosity is intuitively reasonable. 
If the flow and rod radius are regarded as 
fixed, and the porosity is decreased, this 
will give rise to an increase in flow 
velocity and a decrease in passage 
dimensions. Both these effects would be 
expected to contribute to an increase in 
pressure drop, and this is in accord with 
the findings of Figure 6. 

The rod radius is used here as a 
characteristic dimension because it can 
be easily measured and identified in an 
experiment or application. The equivalent 
diameter, = ed/( 1 — «), was also con- 
sidered but not used because it is physi- 
cally somewhat obscure in the present 
instance and did not assist in correlating 
the results. 


Friction-Factor-Reynolds-Number Relationship 

The pressure-drop information may be 
rephrased into a friction-factor-Reynolds- 
number relationship. To begin, the defini- 
tion of the friction factor is 


/ - »/4 < i6 > 

The shear stress may be related to the 
pressure drop, dp/dz. When one refers 
to Figure 5, a force balance yields 

f-r o 0 o dz = — dp- A 2i 
or (17) 

_ _ (dp/dz) A 2 
7*000 

With this Equation (16) becomes 


2 (dp/dz)A >2 

r 0 6 0 pru 2 


(16 a) 
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Fig. 9. Representative velocity contour lines for square array. 


and, after further rearrangement 

- [&(-*)][$] < i6 ‘> 

where N Red is based on the rod diameter 
d(= 2 r 0 ) and is given by 


In writing Equation (166) the authors 
have used the fact that u = Q/A t . 

The first bracket of Equation (166) is 
the pressure-drop-flow ratio, which for a 
given array has been previously shown to 
depend only on the porosity. The second 
bracket is purely geometrical. When one 
uses the definition of e, Equation (15), 
in conjunction with the geometry of 
Figure 5, it is easy to show that for a 
given array the second bracket also is a 
function of porosity alone. 

The friction-factoi^Reynolds-number 
relationship has been computed for each 
array from Equation (166), and the 
results are plotted on Figure 7 as a func- 
tion of porosity. The curve representing 
the triangular array lies higher than that 
for the square array, but intuitive in- 
terpretation is not easy because, in 
addition to the pressure drop, compli- 
cated geometrical factors enter into the 
friction factor. For low porosities it 
would appear that the curves for the two 
arrays are not completely similar in 
shape, but, as previously discussed, 
similarity is not to be expected in the 
low-porosity range. 

For high values of porosity it is 
possible to find a very accurate analytical 
representation of the friction-factor 
results. Under these conditions, where the 
rods are relatively wide apart, the 
velocity distribution around any one rod 
depends very little on the angular 
position. So the velocity can be accurately 
represented by the first two terms of 
Equations (9) and (11). The friction- 
factor-Reynolds-number relationship cor- 
responding to both these abbreviated 
velocity expressions is 

/• A r *„ 

= 8(1 - <*)* (18) 
. 2«* - In «* - (e*)7 2 - 1.5 
where e* = 1 — e. 


By the comparison of the predictions of 
this relation with Figure 7 (which is 
based on the more complete velocity 
solution) it is found that Equation (18) 
is an excellent representation of the 
results for the conditions 

e > 0.80, triangular array 

e > 0.90, square array 

In addition to the curves of the 
present analysis, Figure 7 also includes 
the results of related work. The analysis 
of Emersleben, represented by the 
triangles* and the dashed line, was 
carried out for the square array with 
complex zeta functions used. The nature 
of his solution is such that it can be valid 
only for large values of €, but heretofore 
the actual range of validity has not been 
known. Inspection of Figure 7 shows that 
Emersleben^ results are in close agree- 
ment with those of the current work for 
€ > 0.9, and this may be taken to 
define the limit of validity of his analysis. 

Sullivan's experiments were carried 
out with parallel-oriented filaments of 
goat's wool, blond hair, Chinese hair, 
glass wool, copper wire, and segments of 
drill rod. Each of these materials was 
used separately. The tests utilizing the 
drill rods corresponded to the situation of 
cylinders touching, and the associated 
data points are shown as black circles in 
Figure 7. For the triangular array these 
data are in excellent agreement with the 
present predictions and thereby provide 
support for the theory. For the square 
array the theory was not carried out for 
the situation of cylinders touching 
because of the relatively slower con- 
vergence of the truncated series. However 
inspection of Figure 7 leads to the belief 
. that good agreement would also be 
obtained for this case. For the other tests 
(aside from those with the drill rods) the 
filaments were inserted in a tube in an 
array which was presumably random. 
Hence the experimental conditions corre- 
sponded to neither of the two regular 
arrays studied here. The data, shown as 
open circles in Figure 7, generally fall 
below the analytical curves. It is es- 
pecially interesting that the data con- 
tinue to fall substantially lower than 


♦The points shown on Figure 7 are those reported 
by Carman (£). 



Fig. 10. Local wall shear-stress distribution 
for equilateral triangular array. 



Fig. 11. Local wall shear-stress distribution 
for square array. 


theory even for € > 0.9, where the rela- 
tively large spacing essentially removes 
the effect of rod orientation. This sug- 
gests the possibility of a shortcoming in 
the experimental apparatus, the most 
likely of which would be that the cylin- 
drical filaments were not tightly packed 
against the wall of the bounding tube. 
Such an occurrence would lead to rela- 
tively large, open flow areas near walls 
and a consequent decrease in pressure 
drop and friction factor. 

Velocity Contours 

The distribution of the velocity may 
also be of some interest. The authors 
have confined themselves to representa- 
tive situations, selecting results for a 
spacing ratio s/r 0 of 1.1 to represent 
close packings and those for a ratio 2.0 to 
represent open packings. Dimensionless 
velocity contours (lines of constant 
velocity) are plotted on Figure 8 for the 
triangular array for a typical flow 
element. Inspection of Figure 8a (rela- 
tively large spacing) reveals that the 
velocity contours are essentially circular 
for a sizable region near the rod surface, 
an indication that the neighboring rods 
have little effect there. In Figure 86, 
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where the spacing has decreased, the 
influence of the neighbors extends to 
regions very near the rod surface. Similar 
findings may be observed on the velocity 
contours of Figures 9a and b for the 
square array. 

Variation of Local Wall Shear Stress 

Since the velocity distribution around 
a rod varies with angular position, so 
also will the local shear stress exerted by 
the wall on the fluid. The manner in 
which the shear stress varies is shown 
respectively on Figures 10 and 11 for the 
triangular and square arrays. Using 
symmetry considerations as before one 
need consider only the typical elements 
which are shown cross hatched in Figures 
3 and 4. 

The curves of Figures 10 and 11 
graphically illustrate the role of neighbor- 
ing rods on the flow pattern around a 
given rod. For large spacings, for example 
s/r 0 = 4.0, the local shear stress is 
essentially a constant around the 
periphery of the rod, showing that 
neighbors have little effect on the flow 
pattern. As the spacing decreases, the 
angular dependence of the local-wall 
shear stress increases, testifying to the 
increasing asymmetry of the flow due to 
interference of neighbors. 

As would be expected on physical 
grounds the highest shear stress is 
associated with the location of highest 
velocities (0 = 30 and 45 deg., re- 
spectively), and the smallest shear stress 
is at the location of lowest velocity 
(0=0 deg.). 


CONCLUDING REMARKS 

In the investigation reported here the 
analytical procedure based on truncated 
trigonometric series has been applied to 
symmetrical arrays such as the square 
and the equilateral triangle. However 
the same methods can be applied with 
no essential modification to u asym- 
metrical arrangements, examples of which 
are the rectangular and the isosceles 
triangular with various apex angles. 


APPENDIX 


Integration to Compute Flow Rate 

Equation (13) is the integral for the 
volume flow. For purposes of illustration 
attention is focused on the square array; 
the result for triangular array will be 
given later. 

The expression for the velocity for the 
square array is given by Equation (11). 
The Gj are related to the tabulated con- 
stants S f (Table 1) by 


Gi 


(—dp/dz)s 2 

' 4 jus 41 


(Al) 


For purpose of integration the dimension- 
less variable 


f = r/r 0 (A2) 


is introduced. Then substituting Equa- 
tions (11), (Al), and (A2) into the flow- 
rate integral (13), one gets 



where 0o = 7t/4 has been introduced for 
the square array. Inspection of the con- 
tents of the brace shows that the only 
parameter is the ratio (s/r 0 ), f being 
a dummy integration variable. So the 
assertion of Equation (14) is shown to be 
correct. 

The integrals appearing in Equation 
(A3) can be carried out in a straight- 
forward way, giving 


Qt = (s_Y 

(— dp/d£)r 0 4 Vo/ 



The definite integrals represented by A, 
and /?,• were computed numerically by 
the Kutta-Runge method. 


For the triangular array the authors 
proceed as above, with the exception 
that the velocity is taken from Equation 
(9) with 


Gj 


= A,- 


{ — dp/dz)s 2 
Sjps* 1 


(AS) 


The final result for the flow rate for the 
triangular array is 


qm = M 4 f « 

(- dp/dz)r 0 1 \r 0 / |_ .2 ir \ r 0 



where 


0 , = 
X, = 



cos 6^0 
(cos 0) fl,+3 


dd, 


cos 6j0 
(cos 0) 2-0/ 


de 


NOTATION 

= constants in velocity 
solution 

A\ = solid area normal to flow direc- 
tion 

A 2 = flow area 

d = cylinder diameter, 2r 0 

/ = friction factor, 2f /pu 2 

n = direction of the normal 

N Re = Reynolds number, dup/p 

p = static pressure 

Q = volumetric rate of flow through 

typical element (See Figures 3a 
and 4.) 

r = radial coordinate measured from 

center of cylinder 
r 0 = cylinder radius 

s = half spacing between cylinders 

u = velocity in +2 direction 
u* = reduced velocity defined by 
Equation (3) 
u = mean velocity 

z = longitudinal coordinate 


Greek Symbols 

S, = constants defined by Equation 
(125) 

A,* = constants defined by Equation 

(7 f) 

e = porosity, A%/{Ai + A 2 ) 

€* = solidity, 1 — e 

p = fluid density 

0 = angular coordinate 

0o = angle subtended by flow element 

r = shear stress at cylinder wall 

t = average shear stress 

p — fluid viscosity 
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